Wstęp

Cele

Zbadać częstotliwość oraz skalę opóźnień pociągów w danym okresie oraz zależność ich od sytuacji meteorologicznej. Przedstawić tę zależność za pomocą wizualizacji.

Hipotezy

Przed badaniem robimy przypuszczenie, iż im gorsza sytuacja pogodowa, tzn. im więcej opadów pewnego dnia oraz im mniejsza temperatura, tym opóźnienia będą większe i częstsze.

Zbiór danych

  1. Opóźnienia pociągów

Opis: The dataset contains information about delays on polish railroads published by the Polish State Railways in real time. Data was scraped from https://infopasazer.intercity.pl. It covers the period from the midnight of 2022-05-16 until the midnight of 2022-05-30. Observations were collected every 5 minutes.

Zmienne:

datetime - timestamp at which the sample was collected (Warsaw time)
id - train ID
carrier - the name of the carrier (mainly its PKP Intercity)
date - date of the train departure
connection - beginning and the destination for the train
arrival - planned arrival time at the destination
delay - current estimated delay
name - train station name

Dataset został opublikowany w roku 2022 i do teraz ma nie więcej niż 800 pobrań. Nie znalazłem również żadnych raportów z wykorzystaniem tego zbioru, więc można stwierdzić, że zbiór jest dość unikatowy.

Z drugiej strony, zbiór liczy około 4 mln obserwacji, więc jest dla nas dość reprezentatywny.

  1. Pogoda

Opis: Są to archiwowe dane pogodowe dla okresu, odpowiadającego okresu zbierania danych w pierwszym zbiorze. Za lokalizaję wzięto przybliżone koordynaty geograficznego centrum Polski.

Zarys projektu

W tym projekcie po czyszczeniu zbioru danych oraz konwertowaniu poszczególnych zmiennych do właściwego dla nas typu zajmiemy się badaniem zależności pomiędzy pewnymi typami zmiennych za pomocą np. wyznaczenia korelacji pomiędzy zmiennymi numerycznymi. Znajdziemy też średnie opóźnień dla poszczególnych przewoźników lub relacji, uwzględnimy dane pogodowe, żeby wyłapać tę zależność oraz wesprzymy nasze wnioski odpowiednimi wizualicjami. Pomysłem też jest zrobienie liniowego modelu predykcyjnego, który na podstawie danych pogodowych będzie przewidywał średnie opóźnienie na poszczególnych stacji.

Realizacja

Przygotowanie zbiorów danych

Zacznijmy od importowania zbiorów oraz wszystkich potrzebnych bibliotek

#rm(list = ls())
dftemp = read.csv("delays.csv")
lockBinding("dftemp", globalenv())
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.5
## ✔ ggplot2   3.5.1     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(chron)
## 
## Attaching package: 'chron'
## 
## The following objects are masked from 'package:lubridate':
## 
##     days, hours, minutes, seconds, years
library(ggplot2)
library(fuzzyjoin)
library(corrplot)
## corrplot 0.95 loaded
weather = read.csv("open-meteo-52.55N13.41E38m.csv")

df = dftemp

Zajmijmy się najpierw zbiorem opóźnień

summary(df)
##    datetime              id              carrier              date          
##  Length:3718170     Length:3718170     Length:3718170     Length:3718170    
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##   connection          arrival             delay               name          
##  Length:3718170     Length:3718170     Length:3718170     Length:3718170    
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character
sum(is.na(df))
## [1] 0

Konwertujemy typy danych

df$datetime = strptime(df$datetime, format = "%Y-%m-%d %H:%M")
df$date = strptime(df$date, format = '%Y-%m-%d')
df$arrival = strptime(df$arrival, format = '%H:%M')
split_data = do.call(rbind, strsplit(df$connection, " - "))
df$start = split_data[,1]
df$end = split_data[,2]
rm(split_data)
df = df[, -5]
df$id = as.factor(df$id)
df$carrier = as.factor(df$carrier)
df$delay = as.numeric(gsub(" min", "", df$delay))
df$name = as.factor(df$name)
df$start = as.factor(df$start)
df$end = as.factor(df$end)

Teraz czyscimy zbior danych pogodowych

head(weather)
##           latitude           longitude                elevation
## 1         52.54833           13.407822                     38.0
## 2             time temperature_2m (°C) relative_humidity_2m (%)
## 3 2022-05-16T00:00                12.8                       69
## 4 2022-05-16T01:00                12.0                       70
## 5 2022-05-16T02:00                11.6                       72
## 6 2022-05-16T03:00                10.5                       76
##   utc_offset_seconds               timezone timezone_abbreviation
## 1               7200          Europe/Berlin                  CEST
## 2          rain (mm) surface_pressure (hPa)  wind_speed_10m (m/s)
## 3               0.00                 1016.7                  4.07
## 4               0.00                 1016.7                  3.79
## 5               0.00                 1016.5                  3.04
## 6               0.00                 1016.5                  3.04

Widzimy, że jest niedobrze, pierwsze 2 usunąć

colnames(weather) = c(weather[2, 1], weather[2, 2], weather[2, 3], weather[2, 4], weather[2, 5],weather[2, 6])
weather = weather[c(-1, -2), ]

sum(is.na(weather))
## [1] 0
weather$time = strptime(weather$time, format = '%Y-%m-%d T %H:%M')
weather$`temperature_2m (°C)` = as.numeric(weather$`temperature_2m (°C)`)
weather$`relative_humidity_2m (%)` = as.numeric(weather$`relative_humidity_2m (%)`)
weather$`rain (mm)` = as.numeric(weather$`rain (mm)`)
weather$`surface_pressure (hPa)` = as.numeric(weather$`surface_pressure (hPa)`)
weather$`wind_speed_10m (m/s)` = as.numeric(weather$`wind_speed_10m (m/s)`)

NaNów nie ma, dobrze

Lączymy zbiory

df = fuzzyjoin::difference_left_join(
  df, weather, 
  by = c("datetime" = "time"),     # Adjusting column names here
  max_dist = as.difftime(30, units = "mins")  # Adjust max_dist as needed
) %>%
  group_by(datetime) %>%
  filter(abs(difftime(datetime, time, units = "secs")) == min(abs(difftime(datetime, time, units = "secs")))) %>%
  ungroup() %>%
  select(-time)     # Remove the redundant time column from weather

rm(weather)

Robimy delay factor

df$delay_factor = cut(df$delay,
                       breaks = c(-Inf, 0, 5, 30, 60, Inf),
                       labels = c("0 min", "1-5 min", "6-30 min", "31-60 min", "> 60 min"),
                       right = TRUE)

delay_over_60_by_carrier = df[df$delay_factor == "> 60 min",] %>%
  group_by(carrier)%>%
  summarise("Delays_Over_60_mins" = length(delay_factor))
library(ggridges)

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
delay_summary = df %>%
  group_by(carrier, delay_factor) %>%
  summarise(count = log(n()), .groups = 'drop')

total_counts = delay_summary %>%
  group_by(carrier) %>%
  summarise(total_count = sum(count), .groups = 'drop')

delay_summary = delay_summary %>%
  left_join(total_counts, by = "carrier") %>%
  mutate(count = 100*count / total_count) 

bar_plot = ggplot(delay_summary, aes(x = carrier, y = count, fill = delay_factor)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    title = "Counts of Train Delays by Carrier",
    #x = "Carrier",
    y = "Log(Percentage of delay type)",
    fill = "Delay Factor"
  ) +
  theme_minimal(base_size = 15) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.text.x = element_text(face = "bold"),
    legend.position = "right"
  )

interactive_plot = ggplotly(bar_plot)

interactive_plot = interactive_plot %>%
  layout(
    xaxis = list(tickangle = 45)
  )

interactive_plot
df_by_hour = df%>%
  group_by(format(datetime, "%Y-%m-%d %H"))%>%
  summarize(mean(delay), mean(`rain (mm)`), mean(`temperature_2m (°C)`), mean(`relative_humidity_2m (%)`), mean(`wind_speed_10m (m/s)`))

df_by_hour = df_by_hour[,-1]

correlation = cor(df_by_hour)
corrplot(correlation)